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ABSTRACT 

The shock interaction and evolntion of nova ejecta with a wind from a red 
giant star in a symbiotic binary system are investigated via three-dimensional 
hydrodynamics simulations. We specihcally model the March 2010 outburst of 
the symbiotic recurrent nova V407 Cygni from the quiescent phase to its eruption 
phase. The circumstellar density enhancement due to wind-white dwarf inter¬ 
action is studied in detail. It is found that the density-enhancement efficiency 
depends on the ratio of the orbital speed to the red giant wind speed. Unlike 
another recurrent nova, RS Ophiuchi, we do not observe a strong disk-like den¬ 
sity enhancement, but instead observe an aspherical density distribution with 
~ 20% higher density in the equatorial plane than at the poles. To model the 
2010 outburst, we consider several physical parameters, including the red giant 
mass loss rate, nova eruption energy, and ejecta mass. A detailed study of the 
shock interaction and evolution reveals that the interaction of shocks with the 
red giant wind generates strong Rayleigh-Taylor instabilities. In addition, the 
presence of the companion and circumstellar density enhancement greatly alter 
the shock evolution during the nova phase. The ejecta speed after sweeping out 
most of the circumstellar medium decreases to ~ 100 — 300 km s“^, depending 
on model, which is consistent with the observed extended redward emission in 
[N II] lines in April 2011. 

Subject headings: hydrodynamics — methods: numerical, — novae: general — 
stars:individual: V407 Cygni 
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INTRODUCTION 


The symbiotic recurrent no va (SyRN) V407 Cy g consists of a white dwarf (WD) 


and a Mira-type red giant (RG, 


Munari et al. 


1990l) . In March 2010, the gamma-ray 


e mission from its outburst was detected bv ■Fermz-LAT fLarge Area Telescope. Atwood et 


al. 


2009|), showing the hrst evidence for a gamma-ray nova (lAbdo et ah 


2OI0I ). Another 


ga mma-rav SvRN. V745 Sco. was also detected bv Ferm^-LAT in February 2014 f Cheung et 


al. 


2014J h Since the RGs in symbiotic binaries could produce a high-density circumstellar 


medium (GSM) through their stellar wind, gamma-ray emission in SyRNe could originate 
in high-energy par ticle acceleration that happens when bla st waves pass through the 


high-density GSM flAbdo et al. 


2010 : 


Martin fc Dubus 


2OI3I) 


Re cently, classical novae f CNe) have also been classihed as a new type of gamma-ray 


source (I Ackerman n et ah 


2014J h However, the lack of a high-density GSM in CNe cannot 


explain the origin of gamma-ray emissio n as in SyRNe, requiring a new mechanism for 


gamma-ray CNe. 


Chomiuk et al. 


(120141) studied the high-resolution radio maps of the 


gamma-ray CN V959 Mon, proposing that the fast binary motion in CNe could shape the 
nova ejecta by transferring its angular momentum. Therefore, high-energy particles could 
be accelerated at the interface between equatorial and polar regions. 

Since particle acceleration in CNe and SyRNe is highly associated with shock 
propagation in their nova eruptions, these new discoveries motivate us to study the shock 


evolution in nova systems. Additionally, SvRNe have long been consic 

ered as Type la 

snoernova fSN lai orogenitors (Hachisn & Kato 2001: Kato & Hachisu 

201 

2), although the 

nature of the orogenitor systems of SNe la remains uncertain (Maoz et al. 

2014). However, 

evidence of interaction between SN la ejecta with a GSM has been found in some SNe la 


in late-type spiral galaxies (jPatat et al. 


2007; 


Dildav et al. 


2012 


Silverman et al. 


2013h . 


suggesting that at least some SNe la could originate in symbiotic binaries. 
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The gamma- r ay em ission from the 2010 outburst of V407 Cyg has been studied by 


Martin fc DubusI (120131) through a semi-analytical study of diffusive shock acceleration 


and non-thermal emission in V407 Cyg. They found that the observed gamma-ray light 
curve could be htted by requiring a circumbinary densi ty enhancement (CPE) a round the 
WD. Multi-dimensional hydrodynamics simulations by 


Orlando fc Drake 


20121 ) further 


support the need for a CDE around the WD. They simulated a blast wave passing through 
a dense CSM in V407 Cyg including a sop histicated treatment of thermal conduction and 


radiative cooling fjOrlando 8z Drake 


20121) . Ending that the observed X-ray light curves 


could be reproduced if a CDE exists around the WD. However, they do not include the 
asymmetric effect of binary motion, and the distribution of CDE is artihcially imposed in 
their simulations without modeling the WD-wind interaction in the quiescent phase. 


On the other hand. 


Walder et ah 


(120081 ) performed three-dimensional hydrodynamics 


simulations of another SyRN, RS Oph, from the quiescent phase to a nova eruption. In 
their simulations, they found a disk-like CDE concentrated in the orbital plane, and their 
sh ock evolution is consistent with the observations of the 2006 nova outburst f Sokoloski et 


al. 


2006 


Bode et al. 


20061) . Although s ome similarity between RS Oph and y 407 Cyg has 


been found in their eruptions’ spectra (iMunari et al. 


parameters of the binary systems are quite different (iDrozd et al 


2011 


Shore et al. 


201ll) . the physical 


20131 ). especially for their 


orbital periods (separations). The orbital period in RS Oph is much shorter (Porb ~ 1.25 yr) 
than in V407 Cyg (Porb ~ 43 yr). Therefore, the amount and distribution of density 
enhancement could be different in these two systems. 


In this paper, we investigate the evolution of the SyRN V407 Cyg from the quiescent 
phase to a nova eruption. In the next section, we review and summarize some observational 
facts about the binary system of V407 Cyg and its previous outbursts. The numerical 
methods and the initial setup are presented in Section I3l Our simulation results with 
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different red giant wind loss rates and nova eruption parameters are reported and discussed 
in Section 01 In the final section, we summarize our results and conclude. 


OBSERVATIONS OF V407 CYG 


As mentioned in the above section, the 


a Mira-t ype M6 III RG 


763 day f Kolotilov et ah 


Mnnari et ah 


SyRN V407 Cyg consists of a massive WD and 


19901 ). The Mira-type RG has a pulsation period of 


20031) . From the dust obs curation and t h e long -term optical light 


curve, an orbital period of ~ 43 yr was inferred by 


Mnnari et ah 


(jl990l) . corresponding to 


an orbital separation of 16 AU for a IMq RG and a I.2M0 WD. The distance of V407 Cyg 


Kolotilov et al. 

2003 

1 to 2.7 kpc ( 

Mnnari et al. 

1990) 


1990h. However. Ghomiuk et 


al. (120121) suggest a possible distance > 3.0 kpc using the interstellar Na I D absorption 
lines, since the perio d-luminosity relation i s not well established for the long pulsation 


pe riod of V407 Gvg flWhitelock et al. 


al. 


20121 1 . 


2008h (see a more detailed argument in Ghomiuk et 


V407 Cyg has been monitored in the optical for more than seven decades. Three 


nova eruptions were reported in 1936, 1998, and 2010 ((Schaefer 


2010(1 . although some 


undetected eruptions could have happened in this interval. The March 2010 eruption is the 
most interesting due to the use of modern astronomical instrumentation and the detection 
of gamma-rays. From the measurement of an emission line-width, an ejecta mass 
Mej 10 ^Mq, eru ption energy Ep,^ 10*^^ ergs, and ejec ta speed Vej ~ 3, 200 ± 345 km s ^ 


could be estimated flAbdo et al. 


2010 


Schwarz et al. 


2011 1. An ejecta mass of Mej ~ 10 ^Mq 


could als o reproduce the obser ved velocity evolution, but may underestimate the RG mass 


loss rate ((Ghomiuk et al. 


2012(1 
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In general, the observed mass loss rate M„ind of a Mira giant spans a wide range, 
varying from ^ 10~^Mp yr“^ to ~ 10 ““^Mq yr“^ with an average (M^ind) ~ 6 x 10 “’^Mq yr“^ 


flKnapp &: Morris 


19851) . The radio light curve from the 2010 outburst of V407 Cyg suggests 


a ma ss loss rate ^ 2 x 10 yr ^ ffor v.,a„a = 10 km s ^ and d = 3 kpc: Chomiuk 


et ah 120121 ). but this mass loss rate is one order of magnitude h igher than the esti mated 


value fro m the X-rav light curve ^10 ^ — 10 vr 


& Dubus 


Nelson et al. 


2012 


Martin 


20131) . To obtain better agreement between the X-rav and radio estimates, a 


binary separation wider than 16 AU may be required f) Chomiuk et al 


2012 ). 


The wind speed of a gian t star correlates with the star’s luminosity fjPe Beck et al. 


2010 


Politano fc Taam 


(IMunari et al. 


20 111 ). Taking account of the observed luminosity L lO^Lo 


201lh and the uncertainty of the distance of V407 Cyg, a RG wind speed 
of 10 — 20 km s“^ can be estimated. In addition, a wind speed of ^ 1 0 km s“^ could be 


2010). However, a higher 


also estimated from the optical P Cygni line prohles flAbdo et ahl 
wind speed of 30 km s~^ is estimated fro m the equivalent width of th e Na I D lines in t he 


Nelson et al. 


quiescent phase f|Tatarnikova et al. 


2010 ). 


20031 ) and in the eruption phase (jShore et al. 


2011 


3. NUMERICAL METHODS 

3.1. Simulation Code 


The simula.tion code used in this study is FLASH^I version 4 fiFrvxell et al.ll2000l: Dubev 


et al. 


2008 


Lee 


20131 ). FLASH is a grid-based, parallel, multidimensional hydrodynamics 


code based on block-structured adaptive mesh rehnement ( AMR). The split piecewise 


parabolic method (PPM) solver fjColella fc Woodward 


19841 ) is used to solve for the motion 


^http://flash.uchicago.edu 
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of compressible fluids. We use active particle clouds to represent the WD and the core of 
the RG companion. Particle clouds are utilized in stead of single particles to avoid the fo rce 


anisotropy problem in cloud-in-cell interpolation flRicker fc Taam 


2008 


Panetaf 


2OI2I). 


Gravitational forces bet ween particle clouds and gas and self-gravity are solved using a 


multigrid Poisson solver flRickerl 120081 ). Nuclear burning is ignored in our simulations, since 


it takes place in the envelope of the white dwarf, which is not resolved. We also assume 
that magnetic helds do not affect the shock evolution. 


3.2. Initial Setup 

To model V407 Gyg and its 2010 eruption, we consider a three-dimensional (3D) 
simulation box with a size of 400 AU x 400 AU x 400 AU. Initially, a IMq RG is located at 
the center (origin) of the simulation box, and a 1.2Mq WD is located on the positive x— axis 
with a binary separation of 16 AU, corresponding to an orbital period of about 43 yr. The 
orbital plane is set on the x — y plane and the direction of angular momentum is set to be the 
positive .2—axis. Particle clouds representing the WD and the RG have a radius of 10^^ cm, 
which is slightly larger than three smallest zone spacings (Axmin = 2.93 x 10^^ cm). We use 
9 levels of rehnement based on the magnitudes of the second derivatives of gas density and 
pressure. Each AMR block contains 8^ zones in the 3D simulation box, corresponding to an 
effective uniform resolution of 2048^. To save computation time, we reduce the maximum 
AMR level based on the distance to the center of mass, giving an effective resolution of 
Ax* ~ 0.61 X d at the ith level of AMR, where d is the distance to the center of mass of 
the WD-RG binary. The highest resolution region contains the central 65 AU and the 
hrst forced AMR decrement occurs at d ~ 32 AU. The second AMR decrement occurs at 
d ~ 63 AU, the third at d ~ 127 AU, and so on. Outflow boundary conditions for fluids 
and isolated boundary conditions for the Poisson solver are used. 












We assume the RG has a radius of 2.2 AU, a constant wind speed Uwind = 20 km s“^ 
in the co-rotating frame, and a constant mass loss rate M^md- Therefore, the initial gas 
density p is distributed using Equation [T] 


p{r) = 




wind 




= 1.12 X 10 


-14 




wind 


( '^^wind \ ^ ( ^ —3 111 

(20ten=T) (lAu) 


10-6 Mq yr-1 / V20 km s 

where r is the distance to the RG. We assume a gamma-law equation of state with 7 = 1.1 
and wind temperature T^ind = 7, 000 K to mimic the thermo dynamic behavior of gas 


based on photoionization mo dels of symbiotic binary systems flNussbanmer fc Vogel 


Nussbaumer fc Walder 


19871: 


1993|). 


We perform quiescent-phase simulations up to ~ 2.0 — 2.5 orbital periods, since the 
density distribution becomes quasi-steady in the co-rotating frame after about 1.5 orbital 
periods. During the quiescent phase, we reset the RG wind using Equation [T] for the region 
with r < 5 AU to maintain the RG wind. 


To model a SyRN eruption, we artihcially impose a Sedov-like explosion on the location 
of the WD. The eruption has an ejecta mass Mej, eruption energy Eej, and averaged ejecta 
speed Uej. To avoid grid effects, the eru ption covers a sm all spherical region with a radius 


equal to hfteen smallest zone spacings (jPan et ah 


201211 . Within the erupting region, a 


uniform density distribution and linear velocity prohle in ra dius are used. We assume th e 


kinetic energy is three-fourths of the total eruption energy flDohm-Palmer fc Jones 


199611 . 


The mass and kinetic energy of the GDE within the exploding region and the orbital speed 
of the WD are added to the eruption. Once the SyRN eruption is imposed, we reduce 
the resetting radius for the RG wind to the radius of the RG, and assume a completely 
absorbing surface on the RG. In reality, this assumption is not precisely correct, and a 
reverse shock should be generated when the ejecta reach the surface of the RG. We ignore 
these effects since the purpose of the paper is not to focus on the response of the RG, but 
on the global shock evolution. 
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4. RESULTS AND DISCUSSION 

We perform runs with several different sets of values for the RG wind and SyRN 
eruption parameters, such as Mwmd^-E'ej, and Mej (see Table [T]), and present our simulation 
results in this section. In the hrst subsection, the evolution of wind focused by the WD 
during the quiescent phase is described. The asymmetry due to WD-wind interaction is 
investigated as well. In Section 14.21 we study the shock evolution dependence on different 
RG wind loss rates and different eruption energies and ejecta masses. 


4.1. Quiescent Phase 

We perform two quiescent-phase simulations by varying the mass loss rate from 
^wind = 1 O“®M 0 yr“^ (model M 6 ) to Mwind = yr“^ (model M7; see Table: [I]) to 

study the effects on the GDE. The quiescent phase evolution of both models is qualitatively 
the same, but the gas density distributions scale with their mass loss rates. In this 
subsection, we take model M 6 as the reference simulation and describe the quiescent-phase 
evolution in detail. 

Figure [T] demonstrates a typical evolution of the gas density in the orbital plane during 
the quiescent phase (model M 6 ). Since the initial wind density distribution (see Equation [T]) 
neglects the binary interaction, a turbulent and hot spiral arm quickly forms when the 
RG wind interacts with the WD. This turbulent spiral arm takes about an orbital period 
(~ 40 yr) to reach a quasi-steady large-scale state in the co-rotating frame. We also show 
the corresponding gas temperature distribution in Figure O 

The WD gravitationally deflects matter in its vicinity. Figure |3] shows the mass within 
control surfaces for different radii for our models M 6 and M7. It is clear that the enclosed 
mass increases for the hrst few years and then reaches stability after about 10 yr. Since the 
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RG wind speed is higher than the orbital speed of the WD, the WD experiences a strong 
wind which limits the amount of mass enclosed within the control surfac e. Therefore, w e 


observ e no significant wind-capture disk around the WD as reported by 


Walderetah 


(l2008h . since the wind speed relative to the orbital speed is much higher than that in 


RS Oph. However, we do observe a small, elongated, disk-shaped structure around the WD 
with a size about 5 AU (See Figur e HI). This resu lt is co nsistent with the fast simulation case 


(wwind = 60 km s of RS Oph in IWalder et al.l fl2010l) . since the orbital speed of RS Oph 
('i^orb 15 km s is much faster than that in V407 Cyg (uorb ~ 5 km s ^), giving a similar 
orbital to wind speed ratio of t’orb/'i^wind ~ 1/4. We note, however, that the enclosed mass is 
not gravitationally bound to the WD. As the WD {Rwd ~ 4 x 10® cm) is not resolved in 
our simulations, we cannot determine the mass accreted by the WD. 


For r < 3 — 5 AU, the enclosed mass is enhanced up to ~ 3 — 4 times, independent of 
the radius of the control surface and the mass loss rate. Qualitatively, the models M6 and 
M7 behave similarly, but the density is one order of magnitude higher in model M6 than in 
model M7. In addition, a short-period oscillation (P ~ 2 yr) in the enclosed mass is also 
observed. This prompt oscillation reflects the dynamical timescale around the WD due to 
the interaction between the colliding RG wind and the turbulent flow behind the WD. 


The GDE, including the spiral arm, is concentrated in the equatorial plane. However, 
the GDE does not form a thin disk, but instead shows a mild anisotropy. This is due to 
a relatively large pressure scale height (Pp ~ Cg/VL) in V407 Gyg, where D is the angular 
velocity. Figures [5] and [6] show the gas column density and gas density-weighted temperature 
distribution as viewed from within the orbital plane. The density enhancement is only 
about ~ 20 — 40%, but the temperature can be heated to ~ 20, 000 K. The opening angle of 
the spiral arm is about ~ 30 — 50 degrees. As described above, when the RG wind collides 
with the WD, some gas periodically puffs out from the WD and forms the ring-like density 
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distribution in Figure [5l 

Figure [7] and [H] show the angle-averaged density prohles of models M6 and M7 within 
surface of cones having different inclination angles. Each density prohle is averaged from 180 
uniformly distributed radial rays in a given inclination angle and each radial ray is sampled 
from the surface of the RG to the boundary of the simulation box. From the upper panels, 
the averaged density prohles still follow the 1/r^ distribution, but with about 20 — 40% 
variation for different inclination angles. It is also clear that the density distribution at 
about r < 70 — 80 AU is highly concentrated in the region with 6 < 60°. Beyond this scale, 
the asymmetry is small and only affected at high inclination angles. 

We have also performed a low-resolution run by reducing the maximum AMR level to 8 
and enlarging the hnest zone size to Axmin = 5.86 x 10^^ cm, corresponding to a factor of 2 
larger than the standard resolution. In this low-resolution run, the spiral arm qualitatively 
looks the same as in the standard run, and the averaged density distribution looks the same 
as well, but the turbulence is slightly suppressed due to a higher numerical viscosity. The 
enclosed mass within the same control surface (r < 5Rc) is about 20% higher than the 
standard run. 


It should be noted that in our model we assume a constant RG wind speed during the 
whole quiescent-phase simulation. However, in reality the RG wind is driven by ra diation 


press ure on dust formed above the stellar surface during the Mira gian t’s pulsation fj Willson 


20001 ). The wind acceleration zone could extend to about 5 — IORrg (jBowenlll988l ). which 


is comparable to the binary separation in V407 Gyg. Therefore, the wind speed could be 
less than what we have assumed when it collides with the WD, implying that our estimate 
of the GDE could be low. 
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4.2. Eruption Phase 


After running for two orbital periods in the quiescent-phase simulations, we impose a 
nova eruption on the WD. We perform two eruption simulations with two different eruption 
energies {E = 1.2 x 10^^ erg and E = 1.2 x 10^^ erg) and ejecta masses (Mgj = 10 “’^Mq 
and Mej = IO^^Mq) for each quiescent-phase simulation. Since we assume that the 
averaged ejecta speed is (uej) = 3,000 km/sec and the kinetic ener gy is 

three -fourths of the total eruption energy at the onset of eruption flDohm-Palmer fc Jones 
19961 ). a higher eruption energy implies a larger ejecta mass in that eruption, and vice versa. 
In total, we have four different sets of 3D eruption simulations as described in Table [H In 
Section 14.2.11 we adopt model M6E44 as the reference simulation to describe the overall 
evolution in the eruption phase. We also discuss the shock evolution in all our models in 
Section 14.2.21 


4-2.1. A Qualitative Description of the Evolution after the Nova Eruption 

Figure [9] shows a typical gas density distribution for V407 Cyg in the orbital plane 
after a noya eruption. In this case (model M6E44 in Tabled]), the mass loss rate of the RG 
wind is Mwind = IO^^Mq yr“^, the ejecta mass is Mgj = IO^^Mq, and the eruption energy is 
E = 1.2 X 10^^ erg. At about one week after the eruption, the forward shock sweeps out the 
RG wind and the spiral arm (label A), and forms a reverse shock propagating backward to 
the WD (label B). The two shocks are separated by a contact discontinuity. At ~ 15 days, 
the forward shock reaches the RG. The propagation of the forward shock and the reverse 
shock can be seen in Figure fTOl Note that although we assume the RG completely absorbs 
the forward shock when it reaches the surface of the RG, it does not affect the propagation 
of the reverse shock from the eruption. 
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Subsequently, Rayleigh-Taylor instabilities quickly develop (label C in Figure [9]). 
After the forward shock passing through the RG, the RG wind and the GDE are shocked 
and heated (label D). T hese hot and dense plasma could be the source of X-ray emission 


(lOrlando fc Drakell2012l) . The self-interaction of Rayleigh-Taylor bubbles and the interaction 


of the forward shock with the spiral arm produce additional reverse shocks as well. When 
these reverse shocks collide with each other, hlament-like shocks form at about 70 days 
after the eruption (label E). 


After 70 days, the forward shock has slowed down because of the positive density 
gradient toward the RG, forming a heart-shaped structure. Once the forward shock has 
passed the RG, a high-density tail at the back side of the RG is formed (label F). Within the 
forward shock, several hlaments are formed and destroyed repeatedly due to the interaction 
of the reverse shocks (label G). The turbulent spiral arm is completely destroyed by the 
forward shock. The corresponding gas temperature distribution is shown in Figure [11] 

It is clear that the temperature distribution is highly asymmetric due to the RG wind 
interaction and Rayleigh-Taylor instabilities. High Mach number shocks are associated 
with the forward shock, starting from M > 10 in the hrst week and then decreasing with 
time. Once the Rayleigh-Taylor instabilities have developed, the Mach number decreases to 
M ~ 5 and is maintained at this level for several months. 


A 3D volume rendering of the gas density distribution right before and after 
{t = 137 day) a nova eruption is shown in Figure [T2j The orange color shows the location 
of the forward shock. Several hlaments and Rayleigh-Taylor bubbles can also be seen. 

In Figure [13] and [M] we show the gas column density and gas density-weighted 
temperature distribution in the r — z plane (edge-on view), to demonstrate the aspherical 
evolution of the shock radius due to the GDE. As we described in the previous subsection, 
the GDE makes the density ~ 20% asymmetric at different inclination angles. Therefore, 
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although the ejecta speed is much higher than the orbital speed and the gravitational 
binding energy of the CDE is much less than the eruption energy, the ejecta speed greatly 
decreases in the equatorial direction due to the CDE and the RG wind. The ratio of the 
shock radius in the poleward and equatorial directions is about 1.2 — 1.7, depending on the 
location of the RG. 




Orlando & Drake (120121) have performed simulations of the 2010 eruption of V407 Cyg 


with an artificial CDE. They found that the observed X-ray emission could be reproduced 
if there is a CDE around the WD, but is not required if the outburst energy and ejecta 
mass are near the upper end of the range for these characteristics for classical novae, 
which would be extreme for SyRN. Comparing the results of their best-fitting simulation 
model with CDE, E44.3-NW7-CDE6.3-L40 with our results for model M7E44, we find that 
their forward shock exp ands more spherically and about twice as fast as that found here. 


Orlando fc Drakel (120121) found that most X-ray emission originates from the shocked CSM. 


The observed X-ray e mission reaches a m aximum at around 40 days after the eruption and 


declines after 50 days (IShore et al. 


20 111 ). This corresponds to the time when the RG wind 


and the CDE are heated by the shock in our simulations (see la bel D in Figure 1^. A low er 


shock expansion velocity in our simulations in comparison with 


Orlando fc Drakel (2012) 


may affect the shape of simulated X-ray lightcu rve. However, we note 


averaged ejecta velocity (uej) ~ 9, 000 km s ^ i n 


high based on recent observations (lAbdo et al. 


Orlando fc Drakel (120121) is probably too 


201 


I 


hat the initial 


Since the nova eruption completely destroys the CDE and creates a cavity that 
contains the binary system at the end of the simulation, the subsequent circumstellar 
density distribution will need several years to decades to reestablish the original wind 
profile (Equation [T]). The last recorded nova eruption of V407 Cyg prior to the 2010 
eruption was in 1998. Assuming that there were no undetected eruptions between the 1998 
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and 2010 eruptions and the RG wind has a constant wind speed of 20 km s“^, the newly 
developed wind prohle could only reach to about ~ 50 AU in distance. Therefore, our 
initial wind prohle (Equation [T]) may not be valid at large distances, which may affect our 
shock evolution at late times (t > 3 months for Model M6E44). 

If we reduce the resolution by a factor of 2 as we did in the quiescent-phase 
simulations, we see no notable difference in the shock locations within r < 150 AU, but 
the Rayleigh-Taylor instabilities are less obvious due to a higher numerical viscosity. The 
hlaments produced by the interaction of reverse shocks are slightly less clear than in the 
standard run, but overall the evolution does not change signihcantly. 


4-2.2. Shock Evolutions 

Understanding the shock evolution in nova eruptions is crucial, since the high-energy 
charged particles required for the nonthermal gamma-ray and radio emission are likely 
produced in the shock front. Overall, the simulation results for the shock positions and 
evolution during the eruption phase runs are qualitatively similar, but quantitatively very 
different. Therefore, in this subsection we describe the shock evolution of our four eruption 
simulations in detail. 

Figure 115] to Figure HH] show the shock evolution in all the eruption simulations. At 
the very beginning of a nova eruption {t < 1 week), the eruption is nearly free-expanding 
and the shock location is nearly spherically symmetric. Since the kinetic energy is the 
dominant form of energy at that time and we assume the initial averaged ejecta speed is the 
same in different runs, the shock radius evolution is similar in all cases. However, later on, 
when the ejecta sweep up enough mass, the forward shock slows down due to momentum 
conservation. The decline rate depends on the eruption energy and the RG mass loss rate. 
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For instance, if we lower the eruption energy or the ejecta mass by an order of magnitude 
(case M6E43), the shock location would be slowed down by a factor of ~ 2 at t ~ 200 days 
fFigure [TB]) . Similarly, with the same eruption energy, the shock would propagate much 
faster if the surrounding CDE were less dense fFigure fT7|) . 


Nelson et ah (120121) and Martin fc DubusI (120131) have described a simple semi-analytical 


model for the shock evolution by comparing the circumstellar mass distribution with the 
ejecta mass plus swept mass using momentum conservation. In their model, the density 
distribution is described by a ID wind prohle (Equation H]) plus a CDE. The CDE is 
assumed to be a gaussian distribution which can be described by 


Pcde(’", 6^) — po,CDE exp ( — 


r sin 6 


r cos 6 

I 


( 2 ) 


where Po,cde, I, and b are density and length scale parameters. We perform similar 


calculations by comparing the shoe 
semi-analytical model described by 


location of our 30 hyd rodynamic simulations with the 


Martin fc DnbnsI (2013). 


To hnd the most suitable sets of parameter values, we apply two constraints. First, 
we assume the Po,cde in model M6 is one order of magnitude higher than that in M7, 
i-e. Po,cde(M6)/po,cde(M7) = 10, based on the analysis from Figure [7] and [HI Second, 
we assume the length-scale parameters I and b are equal, since we do not observe 
signiheant asymmetry in the local density enhancement around the WD. By applying 
these two constraints, a best-£t set of parameter values is: for M7, I = b = 20 AU and 
Po,cde = 8 X 10“^® g cm“^; and for M6, I = b = 20 AU and Po,cde = 8 x 10“^^ g cm“^. The 
best-fit shock evolution is shown using the red dotted lines in Figures flSl- [TSl 


The best £t model in 


Orlando fc Drake! (120121 ). which can reproduce the X-ray 


lightcurve, corresponds to their model E44.3-NW7-CDE6.3-L40. It is characterized by an 
ejecta mass of 2 x 10“'^ Mq, eruption energy E = 2 x 10^"^ erg, Po.cde ~ 3 x 10“^® g cm“®. 


and b = I = 40 AU. The CDE in this model is comparable to, in order of magnitude, with 
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our model M7, but is le ss concentrated. However, as described above, the ejecta speed in 


Orlando fc Drakel (120121 ) is much faster than adopted in our work. 


Figure [19] shows the averaged forward shock evolution in the orbital plane and poleward 
regions of all considered models in Table [T] The shock radius evolution can be characterized 
as Tsh oc f", where a is an evolution stage-dependent constant. Depending on the model, a 
is around 0.5 — 0.8 during the first week. At the beginning, the ejecta are nearly spherically 
symmetric, but the angle-averaged forward shock evolves slightly faster in the orbital plane 
than in the poleward direction due to the orbital motion. Later on, the shock interaction 
with the RG reduces the angle-averaged shock speed in the orbital plane, making a decrease 


a little for 10 ^ t ^ 100 days. After about 100 days, the forward shock is 
RG, a nd a increases again. This behavior is consistent with the result in 


ar beyond th e 


Walder et ah 


(120081 ). but for a different time scale due to differences in binary parameters. Furthermore, 


if we increase the RG mass loss rate (from model M7E44 to M6E44 or model M7E43 to 
M6E43) or decrease the eruption energy (model M6E43 to M6E44 or model M7E43 to 
M7E44), a becomes lower. 

The angle-averaged shock speed decreases greatly due to the sweeping out of the RG 
wind. Starting from a value of Uej ~ 3, 000 km s“^ at the onset of the nova eruption, the 
average shock speed in the radial direction drops to v 100 — 300 km s ^ after about a 
year. The evolution of shock speed can be characterized as another power-law relation with 


an index of —| ~ —| (Figure [T9|). The speec 
with the spherically symmetric simulation of 


evolution in Figure [T^ is roughly consistent 


Moore fc BildstenI (1201 21) . In their simulation. 


a better treatment of the radiative cooling during shock expansion is implemented, but only 
in ID. Our shock speed is also comparable with the observed ejecta speed (star symbols in 
Figure [19]) of ~ 2760 km s“^ on day -1-2.3 and 20 0 km s“^ on day -1-196 from the FWHM 


of the broad component oi Ha ([Mnnari et ah 


201ll) . By fitting the simulated ejecta speed 
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with observations, the best-matched models are M6E44 and M7E43. Further more, the 


2011 April observation also shows that the maximum velocity is < 300 km s ^ fj Shore et ah 


201 ll) . which is also consistent with our simulations. The decrement of nova ejecta speed 


with the existence of a CDE in our simulations is also consistent with the observed ejecta 
velocity decrement in RS Oph. 


5. CONCLUSIONS 


We have investigated the symbiotic recurrent nova V407 Cyg from the quiescent phase 
to a nova eruption via three-dimensional hydrodynamical simulations. In quiescent-phase 
simulations, we examined two different mass loss rates of the Mira-type RG. We studied 
the CDE produced due to the wind-WD interaction in V407 Cyg. It is found that the 
induced spiral accretion wake in V407 Cyg is more turbulent and the CDE is less efficient 


than found in RS Oph flWalder et ah 


20081) due to a relatively high wind-to-orbital speed 


rat io. In addition, we o bserve no signihcant wind-captured disk around the WD as reported 


by 


Walder et ah 


(120081) in RS Oph, though the angle-averaged radial density prohle in the 


equatorial plane is about 20% higher than that in the poleward direction. 

The shock evolution in the eruption phase also is studied. It is found that the forward 
shock location is highly dependent on the inclination angle, nova eruption energy, and 
circumstellar density distribution. The forward and reverse shock interactions with the 
RG wind and CDE are also crucial factors in the overall evolution. In addition, the 
shock radius evolution can be characterized by a power law, rsh oc t", where a is about 
0.5 — I.O, depending on model and evolution stage. Our model M6E44 (Mej = I0“®Mq 
and Egj = 1.2 x 10"^“^ erg) and model M7E43 (AUj = I0~’^ Mr7^ and = 1.2 x 10^^ erg) 


show good agreement with the observed ejecta speed from 


Munarietal 


(2011). The shock 


evolution in the presence of CDE in our simulations is similar to what has been observed in 
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V407 Cyg and RS Oph. 

Further work in this series of papers will include investigation of thermal and non- 
thermal emission in V407 Cyg and other symbiotic nova systems at different wavelengths, 
in particular in gamma-ray, X-ray, and radio. Similar analysis can be applied to other RNe 
and CNe as well. 
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Table 1. Simulation Parameters 


Abbreviation M“. . ■ a T’^- a Aft 

Wind wind wind ej ej 

(M 0 yr-l) (AU) (kmsec-l) (K) (erg) (M©) 


Quiescent phase 


Met 

10 -® 

16 

20 

7,000 

— 

— 

M7 

10-7 

16 

20 

7,000 

— 

— 

Eruption phase 

M6E43 

10 “® 

16 

20 

7,000 

1.2 X 10"*® 

10-7 

M6E44t 

10 -® 

16 

20 

7,000 

1.2 X 10'*'^ 

10 -® 

M7E43 

10-7 

16 

20 

7,000 

1.2 X lO'^® 

10-7 

M7E44 

10-7 

16 

20 

7,000 

1.2 X lO'^'^ 

10 -® 


^Reference simulation 


^RG mass loss rate 
^Binary separation 
'^RG wind speed 
^Effective wind temperature 
^Eruption energy 
^Ejecta mass 
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Fig. 1.— Gas density distribution in the orbital plane for model M6 at different labeled 
times. The color scale indicates the logarithm of the gas density in g cm“^. 
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Fig. 2.— Gas temperature distribution in the orbital plane for model M6 at different labeled 
times. The color scale indicates the logarithm of the gas temperature in K. 
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Time (yr) 


Fig. 3.— Mass in the vicinity of the WD companion as determined for spherical control 
surfaces of different radii {Rc = 10^^ cm) as a function of time during the quiescent phase. 
Different line styles indicate different mass loss rates of the RG wind (see Table [T]). 
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Fig. 4.— 3D volume rendering of gas density around the WD for model M6. The white 
region represents the disc-like accretion flow around the WD and the yellow region on the 
left indicates the colliding shock front from the RG wind. The RG cannot be seen in this 
hgure, but the bright region at the bottom shows the location of the RG. 
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Fig. 5.— Gas density projection in the r — z plane perpendicular to the orbital plane for 
model M6 at different labeled times, where the r-axis represents the axis passing through 
the RG and the WD. Note that the origin of the r—axis in this plot is at the center of mass, 
which is different from the simulation coordinate origin. In this hgure, the RG is located on 
the left. The color scale indicates the logarithm of the gas density in g cm“^. 
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Fig. 6.— Similar to Figure [5] but for gas temperature. 
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Fig. 7.— Top: Averaged density profile of model M6 with different inclination angles as a 
fnnction of radins at 82.5 yr. Bottom: The ratio of the density at different inclination angles 
to that along the rotation axis. 
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Fig. 8.— Similar to Figure [7] but for model M7. 
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Fig. 9.— Similar to Figure [T] but for the eruption phase of model M6E44. The labeled times 


indicate the time after a nova eruption. Labels are important features that are described in 
Section 14.21 
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Fig. 10.— Density (solid lines) and pressure (dashed lines) profiles along the line extending 


from the WD to the RG for model M6E44. Different colors represent different times after 
the RN eruption. The red shaded region shows the location of the RG. 
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Fig. 11.— Similar to Figure [2] but for the eruption phase of model M6E44. 
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Fig. 12.— 3D volume rendering of gas density for models M6 and M6E44 at different simula¬ 
tion phases. Left: In the quiescent phase (model M6), the labeled time indicates the simula¬ 
tion time since the start of the simulation. Right; In the eruption phase (model M6E44), the 
labeled time indicates the simulation time after nova eruption. Orange colors in the right 
panel show the location of the nova shock. Red and yellow colors represent high-density 
regions; green and deep blue colors denote low-density regions, (movies of simulations are 
available online). 
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Fig. 13.— Similar to Figure |5] but for the eruption phase of model M6E44. The labeled 
times indicate the time after a nova eruption. 
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Fig. 14.— Similar to Figure [13] but gas temperature. 
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Fig. 15.— Shock radius evolution. Left: in the orbital plane. Right: in the r — z plane 
perpendicular to the orbital plane. Blue contours indicate the shock radius evolution of 
model M6E44 at different times. Black contours show a compariso n of shock radius evolutio n 


with a simple analytical model (denoted by the label MD13) bv Martin fc DnbnsI (120131 ). 


Red dots represent the same simple analytical model, but with the CDE. 
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16.— Similar to Figure [T5] but for model M6E43. 




Fig. 17.— Similar to Figure [T^ but for model M7E44 
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Fig. 18.— Similar to Figure fTSl but for model M7E43 
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Fig. 19.— Time evolution of angle-averaged forward shock radius (left) and velocity (right). 
Different colors represent different eruption simulations in Table [H Solid lines indicate the 
averaged shock radius (velocity) in the orbital plane. Dashed lines show averaged shock 
radiu s (velocity) in the p oleward region. Star symbols represent the observed ejecta speed 


from 


Munarietah 


fcniih . 







